Multigrid Hirsch-Fye quantum Monte Carlo method for dynamical mean-field theory 
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We present a new algorithm which allows for direct numerically exact solutions within dynam- 
ical mean-field theory (DMFT). It is based on the established Hirsch-Fye quantum Monte Carlo 
(IfF-QMC) method. However, the DMFT impurity model is solved not at fixed imaginary-time dis- 
cretization Ar, but for a range of discretization grids; by extrapolation, unbiased Green functions 
are obtained in each DMFT iteration. In contrast to conventional HF-QMC, the multigrid algo- 
rithm converges to the exact DMFT fixed points. It extends the useful range of Ar, is precise and 
reliable even in the immediate vicinity of phase transitions and is more efficient, also in comparison 
to continuous-time methods. Using this algorithm, we show that the spectral weight transfer at the 
Mott transition has been overestimated in a recent density matrix renormalization group study. 

PACS numbers: 71.10.Fd, 71.27.+a, 71.30.+h, 02.70.Ss 



Mott metal-insulator transitions and other effects of 
strong electronic correlations are among the most intrigu- 
ing phenomena in solid state physics They occur 
when the effective electronic bandwidths of cZ or / shell 
electrons become comparable with the local Coulomb in- 
teractions. This most interesting regime is also most 
challenging: here, perturbative expansions (both at 
weak and at strong coupling) and effective independent- 
electron methods such as density functional theory 
within local density approximation (LDA) break down. 

Unfortunately, nonperturbative methods for a direct 
treatment of correlated lattice electrons are either re- 
stricted to one-dimensional (i.e., chain- or ladder- like) 
systems or suffer from severe finite-size and/or sign prob- 
lems. However, a significant reduction of complexity in 
higher-dimensional cases is achieved by the dynamical 
mean-field theory (DMFT) which maps the electronic 
lattice problem onto an Anderson impurity model with 
a self-consistent bath; this mapping becomes exact in 
the limit of infinite lattice coordination Within the 
last 15 years, much insight into strongly correlated sys- 
tems and phenomena has been gained using the DMFT. 
Many of these studies have relied on the Hirsch-Fye quan- 
tum Monte Carlo (HF-QMC) method ^] for solving the 
DMFT impurity problem at finite temperatures. 

The HF-QMC method discretizes the imaginary-time 
path integral into slices of uniform width At and employs 
a Trotter decoupling; this modified impurity problem is 
solved via Monte Carlo sampling of a binary Hubbard- 
Stratonovich field. In principle, arbitrary precision can 
be achieved using HF-QMC in the combined limit of in- 
finitely many updates of the Hubbard-Stratonovich field 
and of vanishing discretization Ar; in this broader sense, 
HF-QMC is numerically exact. However, practical lower 
limits for the discretization exist (primarily due to a scal- 
ing of the numerical effort with (Ar)~^ for fixed tem- 
perature T), so that raw HF-QMC results usually con- 
tain systematic Trotter errors which are much larger than 
the statistical Monte Carlo errors. This implies, in the 



DMFT context, that all observables and phase bound- 
aries depend on the auxiliary parameter Ar chosen in 
the self-consistency cycle. While high accuracy and effi- 
ciency have been demonstrated for a posteriori extrapo- 
lations Ar ^ of certain static observables [^, @] , these 
procedures require special care and experience and may 
fail close to phase transitions. Up to recently, the Trotter 
bias of HF-QMC derived spectral functions could not be 
reduced at all. 

For a long time. Trotter errors could only be avoided 
using fundamentally different impurity solvers that ei- 
ther introduce a logarithmic frequency discretization di- 
rectly (numerical renormalization group, NRG) or repre- 
sent impurity plus bath as a finite cluster or chain [ex- 
act diagonalization (ED), density-matrix renormalization 
group (DMRG)]. Consequently, their results, too, con- 
tain systematic finite-size/discretization errors; in partic- 
ular, these alternatives to QMC yield continuous spectra 
only after numerical broadening. Only very recently, two 
new quantum Monte Carlo impurity solvers have become 
available 0, [1] which are numerically exact in the stricter 
sense that their results are correct within statistical er- 
ror bars, without the need for explicit extrapolations. 
These methods avoid the imaginary-time discretization of 
HF-QMC and rely, instead, on perturbative expansions 
in the interaction and hybridization, respectively, which 
are statistically sampled to arbitrary order. However, 
the continuous- time quantum Monte Carlo (CT-QMC) 
methods are often less efficient than HF-QMC 0. 

In this Letter, we propose a new algorithm for solv- 
ing the DMFT self-consistency equations in quasi con- 
tinuous imaginary time, based on a multigrid implemen- 
tation of the HF-QMC method. This implementation 
removes many of the limitations of HF-QMC; in partic- 
ular, it yields precise and reliable results even at phase 
boundaries. As a first application, we test DMRG pre- 
dictions [9] of the spectral weight transfer at the Mott 
transition - a study that would not have been possible 
using conventional HF-QMC. 
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FIG. 2: (Color online) Scheme: GL free energy in multidi- 
mensional space of hybridization functions {G}- The fixed 
points of conventional HF-QMC (diamonds, squares) are adi- 
abatically connected to the exact fixed point (full circle) only 
for small Ar. In contrast, the multigrid method solves all 
impurity problems at the Ar = fixed point (circles). 



FIG. 1; (Color online) (a) Conventional HF-QMC algorithm: 
Trotter error of impurity Green function G (computed at fixed 
Ar) leads to bias in bath Green function Q, self-energy E, and 
all observables. (b) Extrapolation Ar — » of observable O. 
(c) Multigrid HF-QMC scheme: DMFT iteration with un- 
biased impurity Green function, extrapolated from multiple 
HF-QMC solutions Gati • . ■ Gat„: at self-consistency, G, Q, 
E, and all derived observables are numerically exact. 

Multigrid HF-QMC method - Before we formulate the 
new multigrid approach, let us discuss the flow diagram 
Fig. [H^a) of the conventional HF-QMC method. Starting 
with some guess Sq for the self-energy, the lattice Green 
function G is obtained via the lattice Dyson equation (up- 
per box, here written for the 1-band Hubbard model); 
as a next step, the impurity Dyson equation (middle 
box) is used to determine the bath Green function Q 
which defines the DMFT impurity problem (lower box). 
Conventionally, this is solved using HF-QMC at a fixed 
discretization; the resulting Green function closes the 
DMFT cycle via the impurity Dyson equation. Obvi- 
ously, the Trotter error in the impurity solver affects the 
whole DMFT cycle; thus, at self-consistency, Gat, Qat: 
Eat and all observables deviate from their physical val- 
ues, ideally with corrections in powers of (At)^. Only 
then, numerically exact estimates of observables can be 
extrapolated a posteriori from results of independent HF- 
QMC DMFT runs, as illustrated in Fig.IHb). 

The new multigrid method, visualized in Fig. [Ijc), 
splits the impurity-solving part of the DMFT self- 
consistency cycle into two steps: first, the impurity prob- 
lem is solved ~ at fixed bath Green function Q - for a suit- 
able range of discretizations Ar^ G [Armin, Armax] using 
HF-QMC. In a second step, a numerically exact estimate 
G of the true Green function of the given impurity prob- 
lem is obtained by extrapolation of the multiple result- 



ing Green functions GAri • This highly nontrivial task is 
accomplished using a scheme developed recently (in the 
context of a posteriori extrapolation) [l^l : (i) each of the 
discrete HF-QMC estimates is interpolated with the help 
of a suitably chosen reference model onto a common fine t 
grid; (ii) quadratic (in Ar^ ) least-squares extrapolations 
Ar — > are performed on this grid. In practice, the raw 
Green functions GAri are averaged over 10-20 impurity 
solutions. Thus, for about 10 discretizations Ar,, the 
multigrid method parallelizes to about 100 CPU cores. 
For optimal accuracy, a hierarchy of discretization scales 
and high-frequency cut-offs needs to be maintained [lH; 
details will be presented elsewhere [13] ■ Since the multi- 
grid algorithm closes the DMFT self-consistently with an 
impurity Green function that is (for correctly chosen pa- 
rameters, see below) unbiased, all state variables G, 
and E converge to their numerically exact forms. 

This fundamental advantage of the multigrid HF-QMC 
method is illustrated in the scheme Fig. O Here, the 
solid line depicts the exact Ginzburg-Landau (GL) free 
energy F\U^T,{QY\ in the multi-dimensional space of 
bath Green functions {Q} (for fixed physical parame- 
ters, U^T); its stationary points mark the DMFT fixed 
points [13|]. However, conventional HF-QMC implemen- 
tations converge to the stationary points (squares, di- 
amonds) of generalized GL functionals (dotted lines) 
with additional dependence on Ar. These fixed points 
can vary irregularly or even discontinuously as a func- 
tion of Ar which complicates or prevents extrapolations 
Ar 0. In contrast, the multigrid HF-QMC procedure 
drives the DMFT iteration to the numerically exact fixed 
point: at self-consistency, each finite- Ar HF-QMC eval- 
uation is performed for the bath Green function at which 
F[U,T, {Q}] is stationary (full circle) while the general- 
ized functionals F[U, T, Ar, {G}] are not (empty circles). 
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FIG. 3: (Color online) Double occupancy D — {niUf) at T = 
1/45, U = 5 for metallic (upper set of curves) and insulating 
(lower set of curves) phases: results from conventional HF- 
QMC (diamonds) and from multigrid HF-QMC (circles). 



Benchmark results - For a quantitative discussion, 
let us now consider the one-band Hubbard model with 
semi-elliptic density of states (bandwidth = 4) at 
low temperatures T in the strongly correlated regime. 
Specifically, we will concentrate on the double occupancy 
D — (nin^), an observable which is well-defined for the 
impurity model irrespective of DMFT self-consistency 
and which is best computed at fixed At in both algo- 
rithms. Figure [3] shows results for the coexisting metal- 
lic and insulating phases versus squared discretization at 
T = 1/45 and interaction U = 5. Conventional HF-QMC 
(diamonds) finds an insulating solution only for relatively 
small At < 0.4. While a metallic solution can be stabi- 
lized at arbitrarily large At, the At dependence of D 
is highly nonuniform which restricts the useful range for 
extrapolations At — > (dashed lines), again, to small 
values At < 0.4. In contrast, the multigrid HF-QMC 
raw data (circles) shows regular At dependence even at 
large discretizations which allows for very precise extra- 
polations (solid lines) based on cheap large-AT HF-QMC 
data [iJl- Evidently, the multigrid method is superior. 

However, such extrapolations, as well as all observable 
estimates that can be derived from the Green function 
and self-energy without explicit At ~^ extrapolation, 
are only reliable if the multigrid algorithm has really con- 
verged to the exact DMFT fixed point. Thus, we have to 
check for any bias that might survive from the intrinsic 
parameters, most notably the range [AT,„i,i, AT,„ax] used 
in the internal Green function extrapolation. In the inset 
of Fig. m multigrid HF-QMC results (at [/ = 4, i.e, in 
the metallic phase) are plotted for various At grids; at 
this scale, both the raw data and the fit curves fall on 
top of each other. For a fine-scale analysis, the approxi- 
mate leading Trotter errors have been subtracted in the 
main panel of Fig. (4] Even for this extremely precise raw 
data (with 2 • 10^ sweeps for each At in 10 iterations 
after convergency) , no impact of the grid range can be 
detected for Tmin < 0.25. A slight bias of the order 10~^ 




FIG. 4: (Color online) Inset: multigrid HF-QMC results for 
double occupancy D — {nin-\) at T = 1/45, U = A, using dif- 
ferent ranges [Tmin,Tmax] in the Green function extrapolation. 
Main panel: same data minus leading Ar corrections. 



is visible for the grid range [0.3,0.7]; even the coarsest 
range [0.35, 0.8] of discretizations shifts D only by about 
lO^'*. Similar studies indicate even slightly smaller grid 
range effects in the insulating phase (at U — b) [12| . 
Thus, the multigrid HF-QMC method can be considered 
numerically exact for grid ranges up to [0.3, 0.7] in al- 
most all cases; finer grids or a detailed analysis are only 
needed when extreme precision is sought. This should be 
contrasted with the conventional HF-QMC method for 
which biases in Green functions and self-energies should 
remain detectable down to At w 0.01 - if reasonable 
statistical accuracy could be maintained. 

Spectral weight transfer at the Mott transition - While 
many aspects of the Mott metal-insulator transition are 
well-established by now, including the phase diagram of 
the half-filled frustrated one-band model within DMFT, 
some fundamental questions have remained open. In par- 
ticular, it is still not clear how exactly the spectra change 
across the first-order Mott transition, where the quasi- 
particle peak disappears. Recently, a new scenario has 
been suggested on the basis of dynamical density matrix 
renormalization group (DDMRG) calculations 9]. As 
seen in the inset of Fig. [51 this study finds sharp features 
at the inner edges of the Hubbard bands in the strongly 
correlated metallic phase (solid lines), which are sepa- 
rated from the quasiparticle peak by quite broad pseudo 
gaps; in contrast, the Hubbard bands appear rather fea- 
tureless in the insulating phase (dashed lines) for the 
same interaction U = 5.2, with inner edges that are 
shifted to much smaller frequencies. Obviously this be- 
havior, which has been interpreted in terms of collective 
magnetic excitations in the metal [qI, should be verified. 
One may also ask whether it extends to finite temper- 
atures, i.e., governs the spectral weight transfer at the 
Mott transition, which, for U = 5.2, occurs at T 1/100. 

Estimates of the spectral weight transfer, the difference 
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FIG. 5: (Color online) Inset: DDMRG spectra of coexisting 
metallic and insulating phases at U — 5.2, reproduced from 
[^. Main panel: difference spectra from multigrid HF-QMC 
(broken lines) and DDMRG [3 (solid line). 
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FIG. 6: (Color online) Difference Green functions from QMC 
at finite T (dashed, dotted lines) and extrapolated to T = 
(thick solid line) versus DDMRG data (thin line). Inset: 
extrapolation T — > for t = 1,6, 1.5, 2.5 (top to bottom). 



of metallic and insulating spectra, are shown in the main 
panel of Fig. [51 At this scale, the DDMRG data (sohd 
line) clearly resolves the features discussed above, in par- 
ticular the Hubbard subpeak at |w| « 1.3. The corre- 
sponding multigrid HF-QMC results (broken lines) have 
been obtained, for minimal bias, via Fade analytic con- 
tinuation of Matsubara difference Green functions. They 
are much smoother even at very low T: apart from the 
quasiparticle peak, with a reduced weight, only a shal- 
low negative region at 0.4 < lu < 1.5 is visible, which 
is evidence against the existence of subpeaks. We ob- 
viously cannot exclude that this smoothness is, at least 
partially, an artifact of the ill-conditioned analytic con- 
tinuation; it is also not clear whether the change of shape 
at the lowest temperature T = 1/140, with a closer ap- 
proach of the DDMRG line shape, is significant. In con- 
trast, the QMC difference Green functions shown in Fig. 
in as broken lines (basis for the spectral data of Fig. [5]) 
are precise within linewidth. In the low-r regime that 
we are focusing on, even the extrapolation to T = is 
reliable (see inset of Fig. [6]). However, the transformed 
DDMRG data (thin solid line in main panel) deviates 
markedly from this QMC ground state result (thick solid 
line, with a precision of about linewidth): the DDMRG 
overestimates the differences between metal and insula- 
tor, in this measure, by about 10%, with an even larger 
discrepancy compared to the finite-temperature result at 
T = 1/100, where the metal-insulator transition takes 
place. We conclude that the DDMRG does not capture 
the spectral-weight transfer at the Mott transition quan- 
titatively correctly; possible explanations include the bias 
towards metallic/insulating solutions for odd/even chains 
in the DMRG calculation and, in particular, the deconvo- 
lution technique. In terms of spectra, a large part of the 
DDMRG error can be traced back to its overestimation 
of the quasiparticle peak (since the QMC data of Fig. [5] 



is particularly reliable in this range) ; similar mechanisms 
might have generated the Hubbard-band subpeaks. 

Discussion - We have formulated a new algorithm for 
solving Hubbard- type models within DMFT, based on a 
multigrid implementation of the HF-QMC method with 
internal Green function extrapolation. This technique 
directly yields numerically exact results even at very low 
temperatures, similarly to recently developed continuous- 
time QMC methods. However, our algorithm retains the 
advantages of HF-QMC, i.e., generates Green functions 
with minimal fluctuations, now at about halved temper- 
atures for similar precision and effort. The method ex- 
tends to the general multi-band case, with enormous po- 
tential, e.g., in the context of ab initio LDA-I-DMFT cal- 
culations. While it is too early to judge which of the 
three (quasi) continuous-time methods will ultimately 
prevail, our multigrid method clearly supersedes conven- 
tional HF-QMC for most applications. 

Using this method, we have found inconsistencies in 
the DDMRG scenario of the spectral weight transfer 
at the Mott transition. A definite judgement whether 
the DDMRG results are, at least, qualitatively correct, 
will require additional computing ressources and detailed 
analyses of analytical continuation procedures. 

Stimulating discussions with P.G.J, van Dongen and 
support by the DFG within the Collaborative Research 
Centre SFB/TR 49 are gratefully acknowledged. 
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